This document presents some preliminary analyses characterizing the NDVI of certain places and census tracts throughout the Denver area.

First, load packages used in this Rmarkdown doc.

library(here)
library(tidyverse)
library(sf)
library(terra)
library(raster) 
library(mapview)
library(knitr) #to make sure kable works
library(viridis)
library(scales)
library(lubridate) #for ggplot

NDVI for the Denver area, 2016-2021

We obtained the normalized difference vegetation index (NDVI) from the Landsat-8 satellite at a resolution of 30 meters squared at a roughly 15-day interval between 2016 and 2021 for the following area around Denver:

Define study area

setwd(here("data-processed"))
load("den_metro_bbox_custom.RData")
#around Denver and Jefferson Counties but only as far south as Castle Rock, i.e, not all of Jeff Co to limit data gathered
mv_den_metro_bbox_custom = den_metro_bbox_custom %>% mapview(layer.name = "Study area")
mv_den_metro_bbox_custom

We gathered this data using the rgee package, an R package that facilitates connection to the Google Earth Engine via Python. Please see these two scripts for specific details on the process for gathering and processing the NDVI data in this area:

https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1a_get_landsat_ndvi_denver.R

https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1b_wrangle_check_landsat_ndvi_denver.R

NDVI for one day

Here is the NDVI over the study area on a cloud-free day, July 4, 2021:

setwd(here("data-processed"))
pal_viridis_trunc=viridis::viridis_pal(end=.9) #trunc for truncated
ndvi_den_metro_terr_5_yr = terra::rast("ndvi_den_metro_terr_5_yr.tif")
mv_good_date =ndvi_den_metro_terr_5_yr$`20210704_NDVI` %>% 
  raster::raster() %>%  
  mapview(
    col.regions = pal_viridis_trunc,
    layer.name = "NDVI, July 4, 2021")
mv_good_date

Characterize NDVI of places of interest

First, determine valid (cloud-free) dates

Testing a few places with expected high NDVI

First, pick a few places where we would expect NDVI to be high in the summer on a cloud-free day. If it’s not high, then we can assume the image is bad (i.e., clouds in the way). These test places were determined by examining the cloud-free day (July 4, 2021) above. We chose three areas: an area east of Evergreen, Colorado; a plot in City Park; and a plot in the Indian Tree Golf Club:

(These areas are defined here: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1b_wrangle_check_landsat_ndvi_denver.R)

setwd(here("data-processed"))
load("bbox_evergreen_east.RData")
load("bbox_indian_tree_golf.RData")
load("bbox_city_park.RData")
mv_evergreen_east = bbox_evergreen_east %>%
  mapview(col.regions = "red",  layer.name = "Evergreen East")
mv_indian_tree_golf = bbox_indian_tree_golf %>%
  mapview(col.regions = "red", layer.name = "Indian Tree Golf")
mv_city_park =bbox_city_park %>%
  mapview(col.regions = "red",  layer.name = "City Park plot")
mv_den_metro_bbox_custom = den_metro_bbox_custom %>% 
  mapview(layer.name = "Study area", alpha.regions = .2)
mv_den_metro_bbox_custom+ mv_evergreen_east + mv_indian_tree_golf + mv_city_park

What is NDVI over time for these test places?

Let’s look at the temporal NDVI trend for these three places, without excluding any dates. This is the mean NDVI on that day. That is, each place includes several pixels of size 30 meters squared. This is the average of the NDVI for those pixels for that place on that day.

setwd(here("data-processed"))
load("ndvi_test_places_day_wrangle.RData")
ndvi_test_places_day_wrangle %>% 
  ggplot(
    aes(
      x=date, 
      y=ndvi_mean  
    ))+
  geom_line(aes(colour = test_place_name))+
  geom_point(aes(colour = test_place_name))+
  scale_x_date(labels=date_format("%Y-%m-%d"), date_breaks = "3 months")+
  theme(axis.text.x = element_text(angle = 45, hjust=1))

Hmm, if we look closely, we can see that there are some very low NDVI values for some of these places even in the summer, which is not plausible. Those low values must indicate obstruction by clouds. For example, look at a day when NDVI was measured as very low in City Park despite it being mid-May when we would expect it to be higher. So I’m going to exclude dates like these.

setwd(here("data-processed"))
pal_viridis_trunc=viridis::viridis_pal(end=.9) #trunc for truncated
ndvi_den_metro_terr_5_yr = terra::rast("ndvi_den_metro_terr_5_yr.tif")
mv_bad_date =ndvi_den_metro_terr_5_yr$`20180517_NDVI` %>% #Use 2018-05-17
  raster::raster() %>%  
  mapview(
    col.regions = pal_viridis_trunc,
    layer.name = "NDVI, May 17, 2018")
mv_bad_date

The specific NDVI thresholds for each thresholds are noted in this script (https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1b_wrangle_check_landsat_ndvi_denver.R) and warrant further discussion.

The valid dates are restricted to May, June, July, and August and specifically are:

setwd(here("data-processed"))
load("lookup_date_is_valid_all.RData") #load dataset of valid dates created elsewhere
lookup_date_is_valid_all %>% 
  filter(date_is_valid_all==1) %>% 
  dplyr::select(date) %>% 
  pull()
##  [1] "2016-05-08" "2016-06-01" "2016-06-09" "2016-06-17" "2016-07-03"
##  [6] "2016-08-12" "2017-05-09" "2017-06-02" "2017-06-26" "2017-07-12"
## [11] "2017-07-28" "2017-08-13" "2017-08-21" "2017-08-29" "2018-05-25"
## [16] "2018-06-10" "2018-06-18" "2018-06-26" "2018-07-04" "2018-07-12"
## [21] "2018-08-05" "2018-08-13" "2018-08-21" "2018-08-29" "2019-05-09"
## [26] "2019-06-26" "2019-07-04" "2019-08-21" "2019-08-29" "2020-05-24"
## [31] "2020-06-25" "2020-07-27" "2020-08-12" "2020-08-28" "2021-05-25"
## [36] "2021-06-02" "2021-06-10" "2021-07-04" "2021-07-28" "2021-08-29"

Load places of interest

Load the places-of-interest dataset, which is created here: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/0_read_denver_native_zones.R

setwd(here("data-processed"))
load("places_native_geo.RData")
load("places_native_nogeo.RData")
load("native_places_ndvi_day_nogeo.RData")

Print information about the places of interest, including approximate percent nativity where available.

places_native_nogeo %>% 
  dplyr::select(place_name, native_percent, place_type, place_id) %>% 
  kable()
place_name native_percent place_type place_id
Denver Botanic Gardens, 100% Native 100 native spectrum 1
Plains Conservation Center, 100% Native 100 native spectrum 2
Green Mountain Park, 85% Native 85 native spectrum 3
Hogback along C-470, 75% Native 75 native spectrum 4
Suburban Open Space 1, Chatfield H.S., 50% Native 50 native spectrum 5
Suburban Open Space 2, Columbine Hills Church, 30% Native 30 native spectrum 6
Denver Botanic Gardens at Chatfield, 10% Native 10 native spectrum 7
City Park, 0% Native 0 native spectrum 8
Chatfield Meadow Restoration NA high diversity 9
Chatfield Prairie Garden NA high diversity 10
City Park Greenhouses NA high diversity 11
Denver Botanic Gardens Green Roof NA high diversity 12
Kenderick Lake Xeriscape Garden NA high diversity 13

Map these places by type (nativity spectrum or high-plant diversity).

places_native_geo %>% 
  mapview(zcol = "place_type",
          layer.name = "Place category")

Graph NDVI over time for these places

During valid dates only (above).

For the places with percent nativity values

The first set of polygons sent over.

native_places_ndvi_day_nogeo %>% 
  filter(date_is_valid_all==1) %>% 
  filter(place_type == "native spectrum") %>% 
  ggplot(aes(
    x=date, 
    y=ndvi_50 # median
  ))+
  geom_ribbon( #ribbon around 25th and 75th percentile
    aes(ymin =ndvi_25, ymax = ndvi_75 ),
    alpha=.4
  )+
  ylab("NDVI, Median") +
  xlab("Date") +
  geom_line( size=.7 ) +#note size better than lwd
  geom_point()+
  scale_x_date(breaks = pretty_breaks())+
  scale_y_continuous(
    limits = c(0, NA),  #force axis origin to be zero
    breaks= seq(0,0.8,by = 0.1))+ 
  theme(
    axis.text.x=element_text(angle=60, hjust=1),
        panel.border = element_rect(
          colour = "gray72", size = 0.5, fill=NA))+
  facet_grid(   # facet them
    cols = vars(place_name_fac),
    labeller = label_wrap_gen() #wrap facet labels
  )

For the places with high plant diversity

This is the second set of polygons sent.

native_places_ndvi_day_nogeo %>% 
  filter(date_is_valid_all==1) %>% 
  filter(place_type == "high diversity") %>% 
  ggplot(aes(
    x=date, 
    y=ndvi_50 # median
  ))+
  geom_ribbon( 
    aes(ymin =ndvi_25, ymax = ndvi_75 ),
    alpha=.4
  )+
  ylab("NDVI, Median") +
  xlab("Date") +
  geom_line( size=.7 ) + 
  geom_point()+
  scale_x_date(breaks = pretty_breaks())+
  scale_y_continuous(
    limits = c(0, NA),  
    breaks= seq(0,0.8,by = 0.1))+ 
  theme(
    axis.text.x=element_text(angle=60, hjust=1),
        panel.border = element_rect(
          colour = "gray72", size = 0.5, fill=NA))+
  facet_grid(   
    cols = vars(place_name_fac),
    labeller = label_wrap_gen() #wrap facet labels
  )

Graph NDVI against percent nativity

Among the first set of polygons with values for percent nativity.

native_places_ndvi_day_nogeo %>% 
  filter(date_is_valid_all==1) %>% 
  filter(place_type == "native spectrum") %>% 
  ggplot(aes(
    x=native_percent, 
    y=ndvi_50 # median
  ))+
  geom_point(
    aes(colour=place_name_fac), size = 1.5, 
    alpha=.5 #varying alpha to illustrate density
    ) +
  xlab("Percent native") +
  ylab("NDVI, median") +
  scale_y_continuous(
    limits = c(0, NA),  #force axis origin to be zero
    breaks= seq(0,0.8,by = 0.1))+ 
  scale_color_hue(
    name = "Place name"
  )

Simple summary (median) of NDVI for these places

native_places_ndvi_day_nogeo %>% 
  filter(date_is_valid_all==1) %>% 
  group_by(place_type, place_name_fac) %>% 
  summarise(ndvi_median = median(ndvi_50, na.rm=TRUE) ) %>% 
  ungroup() %>% 
  kable()
place_type place_name_fac ndvi_median
high diversity Chatfield Meadow Restoration 0.4893805
high diversity Chatfield Prairie Garden 0.3686747
high diversity City Park Greenhouses 0.3325997
high diversity Denver Botanic Gardens Green Roof 0.3228098
high diversity Kenderick Lake Xeriscape Garden 0.4737255
native spectrum Denver Botanic Gardens, 100% Native 0.4909220
native spectrum Plains Conservation Center, 100% Native 0.2602929
native spectrum Green Mountain Park, 85% Native 0.3959489
native spectrum Hogback along C-470, 75% Native 0.3483455
native spectrum Suburban Open Space 1, Chatfield H.S., 50% Native 0.4084359
native spectrum Suburban Open Space 2, Columbine Hills Church, 30% Native 0.3175704
native spectrum Denver Botanic Gardens at Chatfield, 10% Native 0.5537543
native spectrum City Park, 0% Native 0.7040144

NDVI of census tracts

The goal of this section is to summarize NDVI for census tracts in Denver County and Jefferson Counties. Before characterizing NDVI of each census tract, we removed bodies of water. NDVI values of water approach -1, and we did not want to “penalize” a census tract for having large amounts of water.

We download census tracts using the tidycensus package in this script: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/0_import_manage_denver_acs.R

We downloaded bodies of water using the osmdata package in this script: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/0_load_denver_osm_parks_water.R

We then removed bodies of water from the census tract here: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/1_merge_tracts_with_water.R

Here is a map of the median NDVI of the census tracts in the area on 2021-06-10. Note that the southernmost portion of the southernmost census tract of Jefferson County is excluded because it was outside of the area from which we gathered data on NDVI.

setwd(here("data-processed"))
load("den_metro_tracts_ndvi_day_wrangle_geo.RData")
load("den_metro_bbox_custom_sf.RData")
pal_viridis_trunc=viridis::viridis_pal(end=.9)  

den_metro_tracts_ndvi_day_wrangle_geo %>% 
  filter(date == "2021-06-10") %>% 
  st_intersection(den_metro_bbox_custom_sf) %>% 
  mapview(
    layer.name = "NDVI, Median",
     col.regions = pal_viridis_trunc,
    zcol = "ndvi_median" 
    )

And the analogous map on 2021-07-04:

Here is a histogram of the distribution of median NDVI for these census tracts on June 10, 2021:

setwd(here("data-processed"))
load("den_metro_tracts_ndvi_day_wrangle_nogeo.RData")
den_metro_tracts_ndvi_day_wrangle_nogeo %>% 
  filter(date == "2021-06-10") %>% 
  ggplot(aes(ndvi_median))+
  geom_histogram()

NDVI of public green spaces

In this section, we examine the NDVI levels of public green space in this area, including city parks and nature reserves. Please see this code for more details on how we gathered the polygons for the greens spaces: https://github.com/michaeldgarber/green-space-denver/blob/main/scripts/0_load_denver_osm_parks_water.R

Here are the parks we downloaded:

setwd(here("data-processed"))
load("den_jeff_co_green_space_public.RData")
den_jeff_co_green_space_public %>% 
  mapview(zcol = "osm_value", layer.name = "Green space type")

Like with census tracts, we removed bodies of water from these green spaces before measuring NDVI. Here is the NDVI on June 10, 2021:

setwd(here("data-processed"))
load("den_metro_green_space_ndvi_day_wrangle_geo.RData")
load("den_jeff_co_geo.RData")
pal_viridis_trunc=viridis::viridis_pal(end=.9) #trunc for truncated
mv_den_jeff_co_geo = den_jeff_co_geo %>%  
  mapview(
    layer.name = "County name",
    color = c("red", "orange"),
    col.regions = c("red", "orange"),
    lwd=2,
    zcol = "county_name_short", alpha.regions = 0)
mv_den_metro_green_space_ndvi_day_geo = den_metro_green_space_ndvi_day_wrangle_geo %>% 
  filter(date == "2021-06-10") %>% 
  mapview(
    layer.name = "NDVI, Median",
    col.regions = pal_viridis_trunc,
    zcol = "ndvi_median" 
  )
mv_den_metro_green_space_ndvi_day_geo+mv_den_jeff_co_geo